library(DEqMS)
library(tidyverse)

df.prot <- read_delim('https://ftp.ebi.ac.uk/pride-archive/2016/06/PXD004163/Yan_miR_Protein_table.flatprottable.txt')

# filter at 1% protein FDR and extract TMT quantifications
df.prot <- df.prot |>
  filter(`miR FASP_q-value` < .01)

# seq(begin, stop, step)
TMT_columns = seq(15,33,2)

df.intn <-  df.prot |>
  select(c(1,all_of(TMT_columns)))

# remove rows with any NA, also work in tibble
df.intn <- df.intn |> na.omit()

df.intn |>
  pivot_longer(-1) |>
  ggplot(aes(name, value)) +
  geom_boxplot() +
  rotate_x_text()

dat.log <- df.intn |>
  mutate(across(2:last_col(), log2))

dat.log |>
  pivot_longer(-1) |>
  ggplot(aes(name, value)) +
  geom_boxplot() +
  rotate_x_text()
  
# if there is only one factor, such as treatment. You can define a vector with
# the treatment group in the same order as samples in the protein table.
cond = as.factor(c("ctrl","miR191","miR372","miR519","ctrl",
                   "miR372","miR519","ctrl","miR191","miR372"))

# The function model.matrix is used to generate the design matrix
design <- model.matrix(~0+cond) # 0 means no intercept for the linear model

colnames(design) <- colnames(design) |>
  str_remove('cond')

# you can define one or multiple contrasts here
x <- c("miR372-ctrl","miR519-ctrl","miR191-ctrl",
       "miR372-miR519","miR372-miR191","miR519-miR191")

contrast <- makeContrasts(contrasts=x,
                          levels=design)
fit1 <- dat.log |>
  column_to_rownames("Protein accession") |>
  lmFit(design) |>
  contrasts.fit(contrast) |>
  eBayes()

# correct bias per protein -------
# based on minimum number of psms
# assign a extra variable `count` to fit3 object, telling how many PSMs are 
# quantifed for each protein
npsm <- df.prot |>
  select(contains('quanted PSMs')) |> 
  as.matrix() |>
  rowMins()

df.prot |>
  select(1, contains('quanted PSMs')) |>
  pivot_longer(-1) |>
  summarise(min = min(value), .by = `Protein accession`) |> pull(min) |> range(na.rm = T)
  
psm.count.table <- data.frame(count = npsm,
                              row.names = df.prot$`Protein accession`)

fit1$genes |> head()
psm.count.table |> glimpse()

fit1$count <- psm.count.table[dat.log$`Protein accession`,"count"]
fit1$count |> head()

fit2 <- spectraCounteBayes(fit1)

# n=30 limits the boxplot to show only proteins quantified by <= 30 PSMs.
fit2 |>
  VarianceBoxplot(n=30,main="TMT10plex dataset PXD004163",xlab="PSM count")

#if you are not sure which coef_col refers to the specific contrast,type
head(fit2$coefficients)

res_mir372_ctl <- fit2 |>
  outputResult(coef_col = 1) |>
  as_tibble()

# volcano plot ---------
# sca means Spectra Count Adjusted
res_mir372_ctl |>
  ggplot(aes(logFC, -log10(adj.P.Val))) +
  geom_point() +
  geom_hline(yintercept = 2) # vignette used p<.001 as threshold, here used fdr<.01
